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A new polarized radiative transfer code Atmospheric Polarization Computations (APC) is 
described. The code is based on separation of the diffuse light field into anisotropic and 
smooth (regular) parts. The anisotropic part is computed analytically. The smooth regular 
part is computed numerically using the discrete ordinates method. Vertical stratification 
of the atmosphere, common types of bidirectional surface reflection and scattering by 
spherical particles or spheroids are included. A particular consideration is given to 
computation of the bidirectional polarization distribution function (BPDF) of the waved 
ocean surface. 

© 2013 Elsevier Ltd. All rights reserved. 


1. Introduction 

Polarization must be taken into account for any exact 
treatment of the problems of light scattering [9, p. 24]. The 
effect of polarization is important both for accurate com- 
putation of the intensity, I, (e.g., GOSAT project, http:// 
www.gosat.nies.go.jp/eng/gosat/page2.htm) and for retriev- 
ing optical properties of the scattering media from the 
measurements of polarization components of the Stokes 
vector if = [/ Q U V] (POLDER-2 instrument onboard 
PARASOL microsatellite http://smsc.cnes.fr/PARASOL). 

Different methods are used for numerical simulations 
of polarized light. Chandrasekhar [9] proposed the solu- 
tion for pure Rayleigh scattering. An adding-doubling 
method was used in the codes developed by Hovenier 
[20], Hansen and Travis [18], and Evans and Stephens [13], 
The method of spherical harmonics was generalized for 
the case of polarization by Benassi et al. [7], Garcia and 
Siewert [14], Ustinov [59], yet this powerful method is not 
widely used today for polarization computations. Other 
solution techniques include the invariant embedding [40], 
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the Monte-Carlo [26,37,42], and successive orders of 
scattering [34] methods. The low order of scattering 
approximations [20,46] are useful for the cases with 
optically thin or strongly absorbing medium. The above 
list of methods is not complete. The reader is referred to 
the monograph by Hovenier et al. [22, paragraph 5.1] for 
more detail. As this monograph is frequently referenced in 
this paper, we will further use an acronym “TPL” (“Transfer 
of polarized light...’’) for convenience. We will also use an 
acronym “ART” (“Analysis of the radiative transfer...’’) for 
another paper [30] frequently referenced herein. In gen- 
eral, the advantages and disadvantages of different meth- 
ods are similar for the scalar and vector cases as discussed 
by Lenoble [33]. 

In this paper, the discrete ordinates method (DOM) is 
used for solution of the transfer problem. DOM was 
developed by Chandrasekhar [9] and references therein. 
It has been used in different vector codes such as 
SCIATRAN [52,53], PStar [47], VLIDORT [56] and others. 
DOM naturally includes the surface reflection in the form 
of Mark [38]. The theoretical basis of DOM as applied to 
the vector radiative transfer equation (VRTE) was devel- 
oped by Siewert [55] and Garcia and Siewert [15]. 

DOM is based on numerical evaluation of the scattering 
integral of VRTE using the Gaussian quadrature. Double 
Gauss rule is often used [57], Given the number of Gaussian 
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nodes (so called “streams”) 2 N, the system of 8 N linear 
differential equations are to be solved (4 N equations per 
hemisphere). By analogy with the scalar case, the VRTE is 
expanded into a Fourier series over azimuthal angle. For 
each of m = 0...M azimuthal harmonics, the system of 8 N 
equations is solved independently. The expansion into 
Fourier series is essentially related to the expansion of the 
phase matrix into Legendre series. In DOM, the number of 
Legendre expansion terms I< and the number of nodes 2 N 
are not independent. The rule 2N=K gives the best result 
according to the properties of the Gauss scheme [58], The 
number K can be very large in case of scattering by large 
particles such as ice crystals, coarse fraction of aerosol or 
hydrosol. This translates into a large number of linear equa- 
tions in the system which increases computational time and 
makes problem ill-conditioned. In the vector case, the number 
of equations is four times as large as in the scalar case and the 
above mentioned problems become more severe. 

The computational time and ill-conditionality can be 
reduced by truncating the number of the Fourier terms (M) 
and the order of Gauss scheme ( N ). Physically it means 


[4,5,28]. The reader is referred to [28] for comparative 
analysis of the MVDOM code and other codes in case of 
coarse aerosol and clouds. APC is a further development of 
MVDOM suitable for simulations and applied analysis: it 
includes vertical stratification of the atmosphere, built-in 
spherical/spheroidal aerosol models [12], and common 
land bidirectional reflectance models or a model of wind- 
ruffled ocean surface as the bottom boundary. More details 
on the new features of APC are provided below. 

The paper is structured as follows: Section 2 defines the 
problem. The anisotropic and regular parts are considered 
in Sections 3 and 4 respectively. The algorithm of compu- 
tation of the surface reflection is described in detail in 
Section 5. Several computational examples are given in 
Section 6. The paper is concluded with the summary. 


2. Formulation of the problem 

The boundary-value problem for the polarized mono- 
chromatic radiation scattered in the plane parallel atmo- 
sphere with no internal sources of energy is given by 


j S (z, p,cp) = jj / q * /_!] P(z,p,(p,p',<p')S(z,p',ip’) dp’ dtp' + Q(z,p, tp,Po, cpoF 

[ S ( 0 ,/<+,<?) = 0 ; S (z 0 ,p_,tp) = l fo" fg R(p,tp,p',tp') S (z 0 ,p',tp')p' dp' dtp' + D(z 0 ,p 0 ,cp 0 ). 


that the radiation field must be a smooth function of the 
view zenith (VZA) and azimuthal angles. In order to 
achieve this goal, the Delta-M method [60] combined with 
the single scattering correction [45] is used. While this 
correction significantly reduces the number of streams N 
and number of harmonics M, it also changes the phase 
function causing substantial error in the aureole area 
[45,31], Recent investigations revealed an error up to 
4.5% [10] caused by truncation of the phase matrix. Thus, 
any vector computations with modified phase matrix 
should be done with caution. 

This paper uses a different approach. Rather than 
changing the phase matrix, it splits the light field into 
the anisotropic and smooth parts. The Stokes vector of the 
diffuse light field is presented in the following form: 

S(z,p,(p) = S A (z,p,tp)+ S R (z,p,tp), (1) 

where S A and S R are the anisotropic [2] and the regular 
parts [4], as defined below, r is the optical depth (OD), p is 
the cosine of the VZA, tp is the relative azimuth. VZA, 9, is 
defined following TPL, Eq. (3.18): p = - cos 9. Eq. (1 ) avoids 
truncation of the phase function. Korkin et al. [30,31] have 
recently investigated approach Eq. (1) for the scalar case 
and showed that it yields a high accuracy solution at low 
orders N and M, similarly to the TMS method [45], but 
with accurate representation of the aureole. 

The approach Eq. (1) is the basis of the new Atmo- 
spheric Polarization Computations (APC) code. In scalar 
case, this approach has been investigated by Romanova 
[51], Irvine [24], and Budak et al. [6], and references 
therein. In vector form, Eq. (1) was used in code MVDOM 


In Eq. (2), a vector with positive cosine p is pointed 
down to the surface, g 0 is the cosine of solar zenith angle 
(SZA), w 0 is the single scattering albedo, z 0 is the total OD 
of the atmosphere, 0<r<r 0 , P(z,p,<p,p' ,cp') is the phase 
matrix, R(p,tp,p',tp') is the bidirectional polarized reflec- 
tance distribution function (BPDF) that accounts for effects 
of polarization. The vector of zeros is noted as 0 . Further 
in the text, the bold symbols indicate square matrices, bold 
symbols with arrows indicate vectors. Q_(z,p,tp,p 0 ,tp 0 ) is 
the free term of the source function 

Q_(z,p,tp,p 0 ,(p 0 )= < ^P(z,p,tp,p 0 ,tp 0 )S 0 exp(-z/p 0 ), (3) 

which represents the direct light beam. For brevity, the 
free term Q_(z,p,tp,p 0 ,tp 0 ) in Eq. (2) (or J in Eq. (30) 
below) is further called "the source function." The term 
D(z 0 ,p 0 ,cpo) describes reflection of the direct solar beam: 

D(z 0 ,p 0 ,cpo) = L/xR(p,tp,p 0 .tpo)Po exp(-z 0 /p 0 ) S 0 . (4) 

Eq. (2) is written in the coordinate system defined by 
Hovenier and van der Mee [21, Eq. (170) and Figs. 2 and 3], 
The solution of Eq. (2) is decomposed into azimuthal Fourier 
series following Siewert [55]. For the case of unpolarized 
collimated beam ( S 0 = [ 1 0 0 0] ) incident on the top 

of atmosphere (TOA), the Fourier decomposition of the 
solution is 


_ M _ 

S(z,p,<p)= £ (2-S 0m )<b(mtp) S ( z,p ), 

m = o 


( 5 ) 
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where 

r f 1 , m = 0 , 

<H<p) = diag[ cos cp cos cp sin <p sin pj; 5 0m = J 

(&) 

The m-th Fourier moment of the source function Eq. (3) 
is given by 

— ^ ^ (On — > 

Q. (T,p,p 0 )= ^-P m (r,p,p 0 )S 0 exp(-r/p 0 ) (7) 

where the m-th moment of the phase matrix 
P(t, p,p', cp, cp') is 


circular CP-basis was proposed by Kuscer and Ribaric [32]. 
The following notation is used here: T c is the matrix of 
transformation from the Stokes basis (SP - Stokes polar- 
ization) to the circular one (CP - circular polarization) and 
T s = T^ 1 gives the inverse transform. Following (Eqs. (1.78) 
and (1.80)) from TPL: 


0 

1 

1 

0 

i 

0 

0 ■ 
1 

: T s = 

■ 0 
1 

1 

0 

1 

0 

o- 

1 

S cp = Tc S sp, 

1 

0 

0 

-1 


-i 

0 

0 

i 

M C p = T c MspT s 

.0 

1 

-i 

0 . 


. 0 

1 

-1 

0. 



( 12 ) 


V rn (x,n,ji')= 2 (2k + l)P^(p)F k (r)P^(p'). 

k = m 


( 8 ) 


Here, P k (p) contains Q™(p), R™(/<) and T™(pi) polynomials, 
related to the generalized spherical functions Pj), n (p) [16,55], I< 
is the maximum considered order of the polynomial. Note, 
that all expansion terms of the phase function are used in 

_ >.m 

Q. 8q. (7), and J (T,p b p 0 ), Eq. (31). 

The scattering matrix and the matrix expansion 
moments (see Eq. (8)), respectively, are written as 
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0 ' 
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0 

, Ffe(r) = 
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0 

0 
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(9) 

where & is the scattering angle. 

Van de Hulst [23, paragraph 5.22], case 6 has formu- 
lated the assumptions behind the symmetry properties of 
the scattering matrix, Eq. (9). The computational expres- 
sions for each element of F^r), Eq. (9), are known from 
[54,41] and TPL, paragraph 2.8. For Rayleigh scattering we 
use Eqs. (3.122) and (3.132) from TPL. 

The scattering matrix Eq. (9) is used for spherical 
particles and spheroids. The APC code uses the model of 
spheres and spheroids described in Dubovik et al. [12]. 

The phase matrix P(z,p,cp,p',cp') takes into account the 
effect of rotation of the reference plane of the Stokes 
vector. The positive angle of rotation is defined following 
TPL, p. 11. Namely, angle of rotation, a, is considered 
positive if rotation is performed in the anti-clockwise 
direction when looking in the direction of propagation. 
Thus, the rotation matrix is (TPL, Eq. (1.51)) 


L(<7) = 


'I 

0 

0 

.0 


0 

COS 2(7 
- sin 2c 
0 


0 o- 

sin 2d 0 
cos 2c 0 
0 1 


L(rr— <t) = L(— cr), 


( 10 ) 


and the phase matrix is defined following TPL, Eqs. (3.7) 
and (3.22), 

P(z,p, Cp,p', Cp') = L(7T-(72)F(r, 0)L( — (7 j ). (IT) 


In order to evaluate the scattering integral, either 
the scattering or the rotation matrix is to be diagonalized. 
The former is possible for the molecular scattering in the 
Chandrasekhar's basis [9, Eq. (162)]. In general case, the 


where M CP and M SP are the representations of an arbitrary 
matrix in CP and SP basis, respectively. 

The definition Eq. (12), used in this paper, is not unique. 
For example, Ustinov [59, Appendix A, Eq. (8)] proposed a 
different equivalent formulation. 

This section provided the main definitions and the 
minimum theoretical background. The next two sections 
will describe solution of Eqs. (1 ) and (2) using the discrete 
ordinates method. 


3. The anisotropic part 


The anisotropic part in vector form has been derived 
and analyzed in detail in Astakhov et al. [2] based on 
properties of the expansion coefficients of S A (z,p, cp) over 
P k m n (p). This section highlights the necessary properties of 
the anisotropic part based on analogy with the scalar case. 
Only a natural (solar) incident beam is considered. 

The anisotropic part, S A (z,p,cp), is computed in the 
coordinate system centered in the solar beam direction 
(see Fig. 1 and ART, Section 2. SAM: evaluation of aniso- 
tropic part). It is convenient to compute S A (z,p,cp) as a 
function of scattering angle 0 = acos (v) where 

v = cos (0) = pp 0 + \/l -p 2 \J 1 -pi cos cp, (13) 

i.e. S A (z,p,cp)= S a (t,v). In the aureole area, i/asl, S A (z,v) 
is close to the accurate numerical solution and almost 
symmetric around the v = 1 direction. Similar to the scalar 
case, it only needs the m = 0 Fourier term. Although the 
S a (t, v) error can be large outside the aureole region, this 
error is a smooth function of angles compensated later in 
the regular part, S R (z,p,cp). 

Following the scalar case [17,33 (Eq. (4.46))], [3], and 
ART, Eq. (3), the VRTE in the small angle approximation is 


Fo- 


dS A (z,p, cp) 
z 


S A (z,p, cp) 



P(T , /(, cp,p\ cp') S A (z,p\ cp’) dp' dcp’ 


+Q_(z,p,cp,p 0 ,cp 0 ), (14) 

where (l(z,p,ip,p 0 ,cp 0 ) is defined in Eq. (3). Eq. (14) is 
subject to the boundary condition S A (0,p,cp)= 0 on the 
TOA only. 

Using the orthogonal properties of the generalized sphe- 
rical functions, the addition theorem (TPL, Appendix B), and 


4 


S.V. Korkin et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 127 (2013) 1-11 



Fig. 1. The right-handed Cartesian coordinate system. OPi is the solar 
beam direction. OP 2 is the direction of observation (scattering). ZOX (the 
same as ZOPi) is the principal plane. PiOP 2 is the scattering plane. ZOP 2 is 
the plane of observation. Positive values of p are measured from the 
positive r-axis (optical depth). The case considered in this figure is 
identical to Fig. 3.3 from TPL. 


the expansion of the VRTE solution [59, Eq. (14)]. 

CP Kmax 2k - i-l „ 

s„ P,v) = Z ^r— Y™ = V)C t (T) (15) 

k = 0 

over the generalized spherical functions (in matrix form) 


Y^“°(/') = diag[ P m, 2 OO p rn,oW 


P m- oW P m,- 2(A) j 


m = Oi 


(16) 


yields the system of equations for the expansion coefficients 
which is similar to that obtained by Ustinov [59, Eq. (18)] 


d C k /dz + C k = Q) 0 X k C k + a ) o exp(-r / p 0 )X k C 0 . (17) 

Recall that Eq. (15) with m = 0 is valid only in the 
coordinate system centered in the solar beam direction. 

In Eq. (17), C 0 is defined by the boundary conditions 
on the TOA 


■Co = T c "S 0 =[0 1/2 1/2 Oj (18) 

and 

Xfc = T C F, ( T S . (19) 

The solution of the VRTE in the small angle approxima- 
tion, Eq. (14), and in CP basis is 

~C k (z) = [exp(-(l-w 0 X k )T//i 0 )-exp(-r/^ 0 )l] "C 0 (t), (20) 


The next equation gives a useful relationship for the 
matrix exponential that follows from T S T C = 1 

T s exp(X k )T c = exp(T s X Jf T c ) = exp(F k ) (22) 

The matrix polynomial is transformed as follows: 

T.sY°(a)T c = P^V) = diag [ Q^V) R° k M Q"m]. 

(23) 

Eqs. (22) and (23) give the expression for the aniso- 
tropic part in the SP basis and in the coordinate system 
centered in the solar beam direction 


S A (r, v) = 2k ^-P° k (p)Z k (T)'S o, 

k - 0 4jI 

where 


(24) 


z k (r) = [exp(-( 1 -ft>oFk )r/ fi 0 )-exp(-z/ n 0 )\\. (25) 

Eq. (25) is computed using the singular value decom- 
position (SVD). Under the natural illumination, Eq. (24) 
contains only two parameters 

~S a (t,p)=[Ia(t,p) Q_a(t:,v) 0 0], (26) 

It follows from Eq. (24) that 1 is defined with 
respect to the plane of scattering PiOP 2 (see Fig. 1). The 
VRTE Eq. (2) is solved with respect to the plane formed by 
the vector of normal, and the view direction (ZOP 2 in 
Fig. 1). In order to match these two coordinate systems, 
rotator Eq. (10) is used exactly as it is used under the 
scattering integral. Following the analogy with the Eq. 
(3.7) from TPL 


S/l(T,i/) = R(ff-(T 2 ) 


z“ —f——P k (v)Z k (T) 

k = o 4 n 


R(-<ti) S o 


= R(-o- 2 )S A (r,i /')=[IaP,p) Q_aP,v) U a (t,p ) 0], 

(27) 


The angle of rotation, <r 2 (Fig. 1 ), is computed using the 
sines and the cosines rule for angles [61]. Eqs. (3.11), (3.12), 
(3.14), and (3.23) from TPL provide the complete set of 
equations for computation of a 2 . Analytically it is conve- 
nient to transform Eq. (15) to the normal coordinate 
system in the CP basis where the rotation matrix has the 
diagonal form, and then to perform the CP-to-SP transfor- 
mation. In the final form, the anisotropic part is expressed 
as follows [5, Eq. (15)]: 

— >• k — >m 

S a (t,/i,^ 0 )= Z ( 2 ~ s 0,m)S A (T,/i,fi 0 ) 
m = 0 

k K m „ 2 /f 4- 1 — 

= Z (2— <5 0 ,m) Z A — ‘K(in f) )P”WZ t (T)Pj , W S o, 

m = 0 k = 0 

(28) 


where 1 is the unit matrix. Eq. (20) is transferred to the 
SP-basis using Eq. (12) 

—).SP —►CP 

S a(t, v)= S a {t,v) = T s S a (t, v) 

K 2k +1 -> 

= T S Z —f— Yt(«/)TcTsC t (r)T c TsCo (21) 

k = 0 4lr 


where <&(mcp) and Z k (r) are defined by Eq. (6) and (25) 
respectively. 

The form of both Eqs. (24) and (28) is valid for the 
vertically stratified case, however Z k (r ) needs to be chan- 
ged. Assuming the case of two layers with z- (1) as the 
optical thickness of the top layer and the optical para- 
meters of two layers defined as ©£, , Fj^ and aij, 2 ’. Fj, 2) 
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respectively, one can write 


matrix form: 


Zfc(r>Ti) = exp(— (1 — <Wq 2) F^ 2) )(t— rt)/^o)CXp(— (1 — gUq ) F^ 1) )ti/^ 0 )— exp(— r/// 0 )l 

(29) 


d S (t) 

dr 


= -M- 1 (l-XW)S(r)+ J (t,mo). 


(32) 


similar to the scalar case [6, Eq. (5.77)]. 

This section provided formulation for the anisotropic 
part using two different coordinate systems and the small- 
angle solution within the forward hemisphere in the 
coordinate system of the incident beam. The anisotropic 
part in the standard coordinate system (with respect to the 
normal) is used next to formulate the bottom boundary 
condition and the source term for the regular part. The 
difference between the total VRTE solution and the 
described anisotropic part is considered in the next 
section. 


4. Regular part 


This section describes the equation and the method of 
numerical solution for the regular part, S R (z,p,<p). Eq. (1) 
is substituted in Eq. (2), and expansion Eq. (5) is used for 
the regular part. There is no particular difference with the 
scalar case at this step (ART, Eq. (9)) except for evaluation 
of the scattering integral in CP basis following Siewert [55], 
The system of 8N ordinary differential equations (i=l... 
2 N, four equations per each value of i ) for the regular part 
is similar to the scalar case (ART, Eq. (16)) 


dS R 
' dr 


= ®q 2 " K 
2 


S R (t,«) 

r 

X (2fc + 1 )WjP^( Mi )F l! (r)P^(pj) S R (t,pj) - 


J (j,Mi,Mo), 
(30) 


where Wj, Mj are the weights and nodes of the Double 
Gauss scheme respectively. Given the number of ordinates 
per sphere 2 N, the number of phase matrix expansion 
terms under the scattering integral is K=2N<K max (see 
Eqs. (15) and (28)). 

Using Eqs. (2) and (14), it is easy to write 

_ >m 

~t m , \ , x dS a (z,mO 

J (t, Mi, f‘o> = (Mo-Mi) ^ 

= E (2k+l)P^(Mi)Z k (r)P^(Mo)~So, 

Mo k = 0 

(31) 

_ >m 

where S A (z,p, tp) is given by Eq. (28). Similarly to (ART, 
Eq. (17)), the recurrence relation for is not used 

here. This simplifies the expression for the source function 
Eq. (31) if compared with [6, Eq. (5.123)]. Eq. (30) is 
independent for each m = 0...M (Eq. (5)). For brevity, the 
upper ( Fourier ) index m is omitted further in the text as 
well as the index “R” that indicates the regular part. Also, 
M~ and m + are used to indicate the directions with negative 
(upward) and positive (downward to the ground) cosines of 
the VZA. The boundary condition on the TOA is (z,p + ) = . 

The boundaiy condition at the surface, S ( t 0 ,m~ ), is discussed 
in Section 5. 

Eq. (30) is solved using SVD and the scaling transfor- 
mation following Karp et al. [25] to improve conditioning 
of the system. For this purpose, Eq. (30) is rewritten in the 


where 


S(T) = [S(T,/i f )] i = 1:2N , J ( r . Mo) — [ J (j,Mi,Mo)/Miii = V2N, 
iij = diag\Mi Mi Mi M 1 , M = diag \ Pi 




P2N 


W = diag[W] ■■■ w 2 n], w, = diag[w,- w, w, w,- 
K 


X = 


E (Zk+))P k (Mi)F k (r)P^(Mj) 


k = m 

Using the notation 

B = M 1 (1-XW) 


J ij= 1 :2JV 


(33) 


and assuming the homogeneous layer, Eq. (32) is solved as 
[25, Eq. (8)] 


— — >• f T ° —> 

-S(0)+exp(Bro)S(ro)= / exp(Br) J 
Jo 


(t) dr. 


(34) 


" l-(0) " 


T- 


R_ 

T_ 

S + (0) 


= 


+ 

T 

R + 


s +(r 0 ) 


J + 


* + 

s -(to) 


Eqs. (20)-(32) from ART give the step-by-step transfor- 
mation of Eq. (32) to the system in the matrix-operator 
form [49]: 


(35) 


where S _(0) and S + (t 0 ) are the solution for the radiation 
reflected from and transmitted through the layer, R + and 
T+ are the elements of the matrix that propagate the 
boundary conditions S + (0) = 0 and S_(t 0 ) through the 
layer. Eqs. (33)-(37) from ART give the step-by-step scaling 
transformation of the source vector J ~j + j . Only two 
features make the processing of the vector case different 
from the scalar one. First of all, the SVD gives complex 
values. They are computed by DGEEV subroutine from the 
LAPACK package. The second feature is that SVD should be 
used not only for Eq. (33), but for Z k {r) in J(t), Eq. (34), as 
well. After that, the integral over r is evaluated analytically. 
The SVD for Z k (r) is performed separately for each 4-by-4 
matrix exponential with a fixed k in the preprocessing. 

In general case, the atmosphere is divided into a 
number of homogeneous layers. The integration in 
Eq. (34) is performed for each layer [25, Eq. (8)]. Using 
two layers, (1) and (2), as an example, one can write 


(36) 


The OD of the upper layer is assumed to be r c , and the 
OD of the second (embedded) layer is assumed to be T 0 -z e . 
Using the continuity condition on the common boundary 
of two layers, one can write 


1 

s_ ( 0 ) 


r ->-(i) i 

J - 


R ( _ 1) 

y(l) " 

r “►(!) i 
S + (0) 

^(1) 
s + (Te) 


_J + J 

+ 

j(l) 

l + 

R^> 

i 

Wl], 

n> 

l 

— ^(2) 
s _ (re) 


r —^(2) 

J_ 


-r(_ 2) 

y(2) " 

r ^ <2) i 
S + (r e ) 

— ^(2) 
s + (to) 


->(2) 

[J + 

H- 

Tf 

R® 

— >(g) 
s _ (to) 


^►d) — ^(2) _>(2) ^(1) 

S + ( z e ) = S + (z e ), S _ ( z e ) = S _ (z e ). 


( 37 ) 
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Eq. (36) is solved with respect to S _ (0) and S + (t 0 ) in 
Eq. (35). After that, the system Eq. (36) can be solved with 
respect to the field inside the media. Currently, the APC 
code computes the reflected and transmitted radiation at 
the interfaces of layers. The height profiles of air pressure 
and temperature are based on MODTRAN [27] and bor- 
rowed from code SHARM [36]. 

The APC code employs the idea of “natural” interpola- 
tion at an arbitrary viewing direction by including the 
required view angles as dummy nodes, p d , into DOM 
scheme with zero weighting coefficients w d = 0 [8[. 

This section briefly described the method of numerical 
computation of the regular part. In general, it follows the 
scalar approach with the exception of complex SVD in the 
vector case. Also, two SVDs are necessary on both sides of 
Eq. (34) in order to perform the integration over OD 
analytically. The use of dummy nodes provides solution 
for the arbitrary directions. The surface boundary condi- 
tion is formulated in the next section. 


5. Surface reflection 


The bottom boundary condition for the regular part is 
derived from Eqs. (1) and (2) 

S r(t 0 ,P-, <p) = - S A (T 0 ,p_,tp) 

r 2 k r \ , 

R(p_,cp,p',tp')[ S A (r 0 ,p’,ip') 


+ - 


i r 1 ’ r 
x Jo Jo 


+ S R (z 0 ,p' ,cp')]p' dp' dtp' 


+ -R(p_,tp,p 0 ,<p 0 )p 0 exp(-r 0 //< 0 )S o- 


(38) 


The matrix-operator technique, Eq. (35), uses the 
azimuthal ( Fourier ) moments, S _(t 0 ), of Eq. (38) 

_>.m — 

S R = - s A (ro,P-) 

+ 


- / R m (ji_,p')[S a (t 0 ,p') + S R (r 0 ,p')]fj' dp' 
xj o 

+R m (/*_,rt)V'o exp(-r 0 /V 0 ) S 0 jn. 


(39) 


The integration in Eq. (39) is performed numerically. In 
the framework of DOM, Eq. (39) is rewritten as 


S-(r 0 )= S R {r 0 ) = — S A (ro)+Rg( S A (r 0 )+ S r {t 0 ))+D b , 
where 


S R (t 0 ) = 


(40) 


- — 

S ) 


- — 

S/iOo.A*) 

— 

S ft (jo,!*-} ) 

; S A (t 0 ) = 

— > m 

S A ( T 0,p}) 

_ ym 


— 

S ft ) 


S a ( T 0? ) 


D B =^ exp(- 


R m Oi>Ao) s o 
R m (Aj,/<o) S 0 

Rm Ojv>Ao) s o 


1 R m (j/j,p+)p+Wi R m (pj,p+)p+w 2 - R m (pf,p+)p+w N 

R B = - 

n ••• ••• 

R m ^,^+)//+w 1 R m (pf,,p+)ft+w 2 ■■■ R m (pf j ,p^)p^w N 

Eq. (40) is substituted in Eq. (35) that now contains 
only two unknowns: S + (t 0 ) and S _(0). The last equation 

in Eq. (35) is solved for S + (r 0 ) analytically. 

In order to perform computations, the Mueller matrix of 
the surface R(p,tp,p',tp') is needed. The following Mueller 
matrix is used to include scalar case in vector computa- 
tions 

R(p,tp,p l ,cp l ) = diag[p(p,tp,p',ip') 0 0 0], (41) 

where p(p,cp,p' ,<p') is the scalar bidirectional reflectance 
distribution function (BRDF). Currently, the APC code 
includes three types of widely used BRDF models for land: 
1 - Lambertian surface; 2a - RPV model [50]; 2b - modi- 
fied RPV (MRPV) model [39] 3 - Ross Thick - Li Sparse 
(RTLS) model [35], The computational subroutines were 
adapted from the code SHARM (see Eqs. (A11)-(A14) and 
(A15)-(A19) from Lyapustin [36]). Eq. (41) is expanded into 
Fourier series 


M 

p(p,p',(p-<p') = (2-8 0im ) X P m (p,p') cos (m(cp-cp')), (42) 

m = 0 

where 

1 f 2 " 

p m (p,p') = 2 n J o p(p,p',tp) cos (mip) dtp. (43) 

Eq. (43) is integrated numerically using the Double 
Gauss scheme. 

Along with the scalar BRDF, code APC includes the BPDF 
model for a wind-ruffled ocean surface (see e.g., [34, 
Section 5[). In this case, the BPDF is 

R (p, p', tp—cp ') = S(p,p\ u, cf>) M(//,/T, cp—cp ')■ (44) 

In Eq. (44) S(p,p',u,c/> ) is the bidirectional shadowing 
function that includes mutual shadowing of waves and 
surface roughness in terms of the probability distribution 
function of the slopes. The roughness of the surface depends 
on the near-surface wind speed, u, for the azimuthally 
independent model by Nakajima and Tanaka [43], The azi- 
muthally dependent Cram-Charlier model [11] includes the 
direction of the wind, rf>. The algorithm of computation of 
the function S(^i,p',u,c/>) is provided in literature (e.g., [36, 
Appendices C and D] for the azimuthally independent and 
azimuthally depended cases respectively). 

The second factor in Eq. (44), M (p,p',ip-<p'), is the 
matrix that includes polarization effects. It contains the 
known Fresnel reflection matrix [47, Eq. (A.l )]. This matrix 
is defined with respect to the plane of incidence and 
reflection of a beam. In general, this plane does not coin- 
cide with the plane defining VRTE solution. In order to 
match these two planes, the rotator is applied to the 
Fresnel matrix. This is similar to the scattering process 
described by the scattering integral in Eq. (2). 

Two methods of computation of M(p,p',i p-ip') are 
known. The first one is described by Eqs. (64)-(79) from 
Mishchenko and Travis [41 ] and uses the amplitude matrix 
that transforms the components of the electro-magnetic 
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field. In this case, the law of reflection of the parallel 
and perpendicular components of the electro-magnetic 
field from a dielectric surface ( Fresnel formulas) must be 
known. Hence, this approach is very convenient for the 
ocean surface. 

Unfortunately, this is generally not the case for most 
natural surfaces including soil, vegetation etc. The Mueller 
matrices of the surface are often measured experimentally 
[48], In these cases, it is necessary to perform the rotation 
of the Mueller matrix explicitly using the rotation matrix 
[1, Eq. (A4)]. This approach does not require the law of 
reflection of the electro-magnetic field, but the angles of 
rotation must be computed explicitly. 

In order to do this, one has to take into account that the 
Mueller matrix of the surface, including the rotators, is 
written with respect to the surface normal [41 ] which is 
the opposite to the direction of the positive cosines of the 
VRTE, Eq. (2). For the surface reflection, both incident and 
reflected angles are often assumed to be within [0, 90) 
degrees [50] which is also different from the definition of 
angles used in the VRTE. 

In order to avoid the above mentioned mismatches, the 
explicit relation for M(p, pi ' , cp-q>') is provided in this paper. 
He et al. [19] have provided the reflection-transmission 
Mueller matrices of the facet (f 2.4.1 and 2.4.2 of their 
paper), but they did not explicitly discuss the problem 
of Fourier expansion of the Mueller matrix. The process 
of reflection is considered similarly to the process of 
scattering with the obvious exception, that for reflection 
the directions of incidence and “scattering” always belong 
to the same hemisphere. In order to match the coordinate 
systems of the VRTE and the reflecting surface, we use one 
normal direction as defined in Section 1. It means that the 
signs of the angles of rotation during reflection are equal 
to those during scattering (see Eq. (11)) 

M(/f, p', cp-rp') = L(7r-t7 2 )Fj{(0)L(-ffi ), (45) 

where 8 is a half of an angle complementary to the 
“scattering” angle Eq. (13). The elements of the Fresnel 
matrix 


Eq. (49) allows to compute the elements of Fr(i 9) in 
Eq. (46). The angles in Eq. (47) are computed using 
Eqs. (3.11), (3.12) and (3.23) from TPL without any mod- 
ifications. Here the particular attention is paid only to 
Eq. (3.23) from TPL 

sintri _ siner 2 _ sin(y-y') 

\J 1 -p 2 \J1-p' 2 -/l- cos 2 © 

where L’Hopital's rule is used in Eq. (49) for \p\ = 1 or 
\u'\=\. Finally, no rotation is performed for <p-cp' = 0, n 

(TPL, p. 71). 

The azimuthal expansion of Eq. (48) is similar to the 
azimuthal expansion of the matrix Eq. (11). Using Eq. (18) 
from Siewert [55] and the definition 


C m = cos (mtp); S m = sin(mr/)), 



(51) 

one can write 

r /\/f m 

r 

C A/f m 

0 ' 

M 

M(p,p',ip-cp')= £ Q.-S 0m ) 
m = 0 

r 

<: 

r M m 
'-m‘W22 

J m lvI 32 

C A/T m 

r 

y -m lvI 33 

0 

0 


0 

0 

0 

r 

K-m IVI 44 


M 

= it (2-<5o, m )I ip-cp'). (52) 

m = 0 

Eqs. (48) and (52) provides the law of Fourier expansion 
of the elements of the Mueller matrix. For example, M™ is 
computed using Eq. (42). Lenoble et al. [34, Sections 5.2.1 
and 5.2.2] provided detailed analysis of Eqs. (44)-(48). 
Similar to the phase matrix, they used analytical integra- 
tion over azimuth to get the Fourier expansion moments 
M m (pi,fi',cp-ip'). This approach is more complicated than 
the numeric integration over azimuth and still requires 
numeric integration over zenith angle to compute expan- 
sion moments ([34, Eq. (45)]). 

To conclude this section, we note that APC uses the 
postprocessing correction of the direct beam reflected 
from the surface and transmitted to the TOA. The correc- 
tion term is expressed as 


Fr(6>) = 


r+ 

r_ 

0 

0 


r_ 

r+ 

0 

0 


0 0 " 

0 0 
r 0 0 ’ 
0 r 0 


(46) 


are known [47, Eqs. (A1)-(A3)[. Using the notation 
c 1-2 = cos(2<t 1i2 ); s 1-2 = sin(2(r u ), (47) 


y p JVJ 

At<m(P-,<P) = — £ ( 2 -<5o,m )R m (P_ , Po . ‘Po ) 

n L m = 0 J 

xexp( ^'° (53) 

V PoP- J 

where R m (p_,<p,p 0 ,tp 0 ) is defined by Eq. (52). The correc- 
tion term, A T0A , is applied to the resulting Stokes vector at 
the TOA 


and Eq. (45), one can write 

r + r_Cj -r_S] 0 " 

r_c 2 r + cic 2 -roSis 2 -r + sic 2 -r 0 cis 2 0 

M (u, u' = _ 

r_s 2 r + c,s 2 + r 0 s 1 c 2 -r + s^s 2 + ^ 0^2 0 

0 0 0 r 0 

(48) 

Further on, the following definition {p < 0, <p ) and 
{p' > 0, <p ' ) is used for the direction of reflection (observa- 
tion) and incidence, respectively. The angle 28 is 

cos 28= cos ( ic-9 ) = —pp'~ \/ 1 -p 2 a/ 1 -p' 2 cos 1 p. (49) 


s TOA(r = 0,p_,tp)= S (t= 0,p_,(p) 4- A toa(F-,<P) (54) 

and does not affect transmission of radiation through the 
atmosphere. Eq. (53) exactly describes the surface signal, 
transmitted through the atmosphere without scattering. 
The forward peak is included in the small-angle part and 
does not require correction. 

This section provided the detailed formulation for the 
surface reflection. The algorithm of computation of the 
Mueller matrix that includes the explicit rotation of the 
reference plane is described in the coordinate system of 
the VRTE Eq. (2). All necessary expressions are provided 
either explicitly or by means of the direct reference to the 
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Fig. 2. The relative deviation, %, between APC and SCIATRAN for / (left column), Q (middle column), and U (right column) components of the Stokes vector. 
Solar zenith angle 60 , relative azimuth 0 (red lines) and 45' (blue lines). View zenith angles from 0“ to 80 , with 10 step (reflected radiation). Rayleigh 
atmosphere with optical thickness 0.1 (top row) and 0.5(bottom tow). Single scattering albedo 0.9, no depolarization. Ruffled ocean surface was simulated 
using Nakajima-Tanaka model. Wind speed 1 m/s at 10 m above the surface. (For interpretation of the references to color in this figure legend, the reader is 
referred to the web version of this article.) 


published equations. This concludes the description of the 
theoretical background of the APC code. In the last Section 
6 some numerical results are provided. 

6. Numerical results 

The code APC, written in FORTRAN, was thoroughly 
tested against its Matlab prototype, MVDOM, which in turn 
was validated against SCIATRAN, Pstar, and other codes 
involved in the benchmark test described in Kokhanovsky 
et al. [28]. For a homogeneous and vertically stratified 
media and for black and Lambertian reflecting surfaces the 
code was tested in Korkin et al. [29] showing good (five 
significant digits) agreement with RT3 [13]. Code SHARM 
[36] was used to validate the scalar surface BRDF models 
adapted in code APC. 

To test the ocean BPDF model given by Eqs. (43) and 
(48), and atmosphere-ocean interaction, we made a direct 
comparison with code SCIATRAN [53]. The SCIATRAN soft- 
ware package along with a detailed User's Guide is avail- 
able via the webpage of the Institute of Environmental 
Physics (IUP), University of Bremen: http://www.iup.phy 
sik.uni-bremen.de/sciatran. As an example, we present 
computations of the reflected radiation over the ocean in 
Rayleigh case for the geometry 0 o =6O", relative azimuths 
tp = 0 (dash lines) and 45° (solid lines), and the VZA from 


0 to 80 with 10 degrees step. The following atmospheric 
parameters (OD t=0.1, and 0.5, SSA m o =0.9, and zero 
depolarization factor) and the wind speed 1 m/s with 
Nakajima and Tanaka [44] model were used. The relative 
difference between the APC and SCIATRAN results for I, Q, 
and U (for nonzero azimuth only) is shown in Fig. 2. The 
relative deviation between the two codes was found to be 
negligible. 

To illustrate functionality of the developed code, the 
next example shows the sensitivity of the intensity, 1, and 
the degree of polarization 

P =s Jq 2 + U 2 + V 2 /I (55) 

to the shape of particles represented by spheres and 
spheroids. The simulations were performed for an idea- 
lized two-layer atmosphere illuminated at SZA=45 D . The 
top layer is assumed to be purely Rayleigh with optical 
depth t r =0.3, and the bottom layer consisting of the non- 
absorbing aerosol only. The surface bidirectional reflec- 
tance was simulated using RPV model with parameters 
typical of savannah [36]. The aerosols are represented by 
spherical particles and spheroids with integrated aspect 
ratio close to 3. Fig. 3 shows the phase matrices for spheres 
and spheroids, respectively. Note, that all diagonal 
elements, a^ 34 , of the scattering matrix contain scattering 
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Scattering angle, degrees 



Fig. 3. The elements of the phase matrices for spheres (red line) and spheroids (blue line) used for computation. See details in Section 6. (For interpretation 
of the references to color in this figure legend, the reader is referred to the web version of this article.) 


anisotropy. Computations were performed for two values of 
aerosol optical depth, t a =OA, 0.5 and the geometry of 
VZA=180°, 150', 120° (transmitted radiation), and 0°, 30°, 
60° (reflected radiation) for the full range of azimuthal angles 
with 5" step. The blue, green and red colors represent different 
VZA (e.g., to VZA=0°, 30 , 60 for the case of reflection). 

The results for the total intensity, /, and for the degree 
of polarization, P (see Eq. (55)), are shown in Fig. 4 for both 
reflected and transmitted radiation. The solid and dashed 
lines represent spherical aerosols and spheroids, respec- 
tively. In the considered case, the figures clearly show that 
polarization is more sensitive to the shape of the particles 
even if depolarizing surface is assumed. In case of nadir 
observation, sensitivity of the total intensity to the parti- 
cles' shape is negligible (Fig. 4a and e) whereas sensitivity 
of polarization remains considerable (Fig. 4b and f). The 
same holds true for the transmitted radiation. 

7. Conclusion 

In this paper, we presented a new code, APC, that uses 
the decomposition of the diffuse light field into a regular 
and anisotropic parts. Contraiy to approximations that 
truncate the phase matrix, this approach does not modify 
the scattering law. As in the scalar case, the regular part of 
the signal, which requires the numerical solution, is 
smoothed in both zenith and azimuth angles after sub- 
traction of the anisotropic component. The anisotropic 
part is computed analytically, while the discrete ordinates 
method is used for numerical computation of the regular 
part. The code is developed in Fortran 90/95 and uses 


LAPACK package to perform the singular value decomposi- 
tion and matrix operations. 

For the users' convenience, the APC code has built-in 
atmospheric profiles and several common models of the 
surface bidirectional reflectance. This paper paid a parti- 
cular attention to the problem of reflection of polarized 
light from the waved ocean surface. We have provided the 
explicit formulas and direct references to the equations 
that describe the reflection of polarized light from any 
surface, for which the Mueller matrix is known, while the 
law of reflection of the components of the electro-mag- 
netic field may be unknown. Although the equations are 
not new, we are not aware of the detailed description of 
the complete model provided here. In our approach, an 
attempt was made to keep the same coordinate system for 
the VRTE problem and the boundary condition. The code is 
available at ftp://climatel.gsfc.nasa.gov/skorkin. The kernel 
databases for fast spheric and spheroidal scattering com- 
putations are also provided by courtesy of Dubovik [12]. 
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Reflection, x = 0.5 




Transmission, x = 0.5 




Azimuth angle, degrees 


Fig. 4. The values of the total intensity, /, and the degree of polarization, P, in case of spheres (solid lines) and spheroids (dots) as functions of the relative 
azimuth. VZA=180° (red), 150° (green), 120° (blue) for the transmitted light (c, d, g, h) and VZA=0° (red), 30° (green), 60° (blue) for the reflected light (a, b, 
e, f). The results are shown for two different optical depths of the aerosol layer. (For interpretation of the references to color in this figure legend, the reader 
is referred to the web version of this article.) 
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